Negative radiation pressure in metamaterials explained by light-driven atomic mass density rarefication waves

The momentum and radiation pressure of light in negative-index metamaterials (NIMs) are commonly expected to reverse their direction from what is observed for normal materials. The negative refraction and inverse Doppler effect of light in NIMs have been experimentally observed, but the equally surprising phenomenon, the negative radiation pressure of light, still lacks experimental verification. We show by simulating the exact position- and time-dependent field-material dynamics in NIMs that the momentum and radiation pressure of light in NIMs can be either positive or negative depending on their subwavelength structure. In NIMs exhibiting negative radiation pressure, the negative total momentum of light is caused by the sum of the positive momentum of the electromagnetic field and the negative momentum of the material. The negative momentum of the material results from the optical force density, which drives atoms backward and reduces the local density of atoms at the site of the light field. In contrast to earlier works, light in NIMs exhibiting negative radiation pressure has both negative total momentum and energy. For the experimental discovery of the negative radiation pressure, one must carefully design the NIM structure and record the joint total pressure of the field and material momentum components.

Negative-index metamaterials (NIMs) exhibit a number of extraordinary optical phenomena ranging from the inverse Doppler effect to beating the diffraction limit with superlenses [1][2][3][4][5][6][7][8][9] . The possibility of NIMs was originally hypothesized by Veselago in 1960s 10 . Among the astonishing phenomena predicted by Veselago, the reversal of the momentum density of light is still experimentally unverified. The reversal of the momentum density of light to be opposite to the electromagnetic (EM) energy flux density predicted by Veselago is counter-intuitive. It is well known that any field, particle, or quasiparticle, having positive energy propagating in one direction, must have its momentum in the same direction. Accordingly, the particle must fulfill the standard momentum per energy ratio p/E = v g /c 2 of relativistic mechanics 11 , where v g is the group velocity and c is the speed of light in vacuum. Therefore, if the total momentum of the coupled field-material state of light in NIMs is truly negative, the light in NIMs should be described by negative-energy quasiparticles, illustrated in Fig. 1. The energy of these quasiparticles would not be equal to the known photon energy ω 0 , which is positive. This raises a question how these quasiparticles emerge when a photon is transferred from vacuum to a NIM in such a way that both the photon energy and momentum are transferred to the quasiparticle and the interface layer of the NIM crystal. The present work answers these questions by providing a detailed position-and time-dependent theory and the related computer simulations of the energy and momentum of light in NIMs.
In recent literature on optical forces in inhomogeneous materials like NIMs, it is typical to assume the Maxwell or Helmholtz stress tensor and further take the time average of the related optical force density [12][13][14][15][16][17] . Homogenization of material parameters and fields may also be applied 12,13,[18][19][20] . On the experimental side, the optical forces on metasurfaces have been very recently applied to drive microvehicles 21 , but their use has not yet reached phenomena of light propagating inside NIM structures. The non-averaged exact position-and timedependent optical force and momentum densities of light in NIMs can be obtained by solving simultaneously Maxwell's equations for the field and the Newtonian dynamics of the solid, neutral material. The full picture of the energy and momentum fluxes associated with optical fields in NIMs can be obtained after solving the coupled dynamical equations. 1 Photonics Group, Department of Electronics and Nanoengineering, Aalto University, P.O. Box 13500, 00076 Aalto, Finland. 2  www.nature.com/scientificreports/ In this work, we use the covariant mass-polariton (MP) theory of light [22][23][24][25][26][27] to analyze the energy and momentum of light in NIMs. We also discuss experimental discovery of reduced, p < p 0 , or even negative, p < 0 , radiation pressure, where p and p 0 denote the components of total momentum of light in the NIM and in vacuum in the direction of propagation. The MP theory describes fields by using conventional Maxwell's equations, which explain the experimentally verified negative refraction and the inverse Doppler effect. However, our approach is fundamentally different from previous theories of negative radiation pressure in NIMs since it describes light in NIMs as coupled field-material MP states of Fig. 1. In the MP states of NIMs, the positive EM energy and momentum densities are accompanied by positive or negative field-driven energy and momentum densities of the material constituting an atomic mass density wave (MDW). We also show that, in the quantum picture, light in NIMs must be described as MP quasiparticles, whose total energy and momentum can become negative. The theory applies to both hypothetical homogeneous and conventional inhomogeneous NIMs structures. It shows that the total momentum of light in NIMs is typically smaller than the momentum of light in vacuum, and depending on the structure, it can even become negative. This contrasts to normal materials, for which the total momentum is always greater than its vacuum value. Our results demonstrate that the MP theory provides a versatile tool for engineering of optical forces in arbitrary NIMs structures and for designing of measurements, where the reduced or negative radiation pressure can be experimentally verified.

Results
Momentum and energy densities of the EM field. We approach the description of light in arbitrary lossless photonic crystals from the classical field theory point of view assuming that the material is at rest in the absence of light. The energy flux density described by the Poynting vector S = E × H , where E and H are the electric and magnetic fields, is known to be continuous over material interfaces. This leads to a continuous EM contribution to the momentum density of light, given by 28,29 As discussed in "Methods", the energy density of the EM field in the MP theory of light in dispersive media, solved from the conservation law of energy, is given by 23 Here ω 0 is the central angular frequency of the narrow frequency band, and ε = ε(ω 0 ) and µ = µ(ω 0 ) are the permittivity and permeability of the material at ω 0 . For a narrow frequency band, the permittivity can be approximated by the first two terms of its Taylor series ε(ω) ≈ ε(ω 0 ) + ∂ε(ω 0 ) ∂ω 0 (ω − ω 0 ) . A similar expansion is valid for the permeability. The angle brackets in Eq. (2) are defined by The brackets here and below are also equal to the time average of the pertinent quantities over the harmonic cycle, e.g., �E 2 (r, t)� = 1 T T/2 −T/2 E 2 (r, t + t ′ )dt ′ , where T = 2π/ω 0 is the length of the harmonic cycle.
Optical force density. The electric and magnetic fields E and H and corresponding flux densities D and B of light in NIMs can be solved from Maxwell's equations. What is less well known is the dynamics of atoms in NIMs under the influence of the instantaneous optical force density. The microscopic form of the optical force density in the MP theory of light accounting for the gradients of material parameters in inhomogeneous dispersive materials is given by 23 (1) (2) The optical field is incident from the left. Inside the NIM crystal, it forms a coupled quasiparticle state with atomic displacements driven by the optical force density. When the field enters or exits the NIM crystal, the interfaces of the crystal take recoil momenta due to the momentum change of the EM wave in vacuum and the coupled state in the NIM crystal. In contrast to normal materials, the energy of the coupled state is negative, and its momentum is directed backward. The directions of the interface momenta are also reversed from what is observed for normal materials. Here I is the 3 × 3 unit matrix and ⊗ denotes the outer product of vectors.

Momentum and energy densities of the atomic MDW and the coupled MP state of light.
When light is interacting with a lossless material, the dynamics of the material is described solely by the Newton's equation. Assuming that there are no other forces except the optical force, including both interface and bulk forces, Newton's equation of motion for the material is given by 22,30 Here p a is the momentum of a single atom and n a is the number density of atoms. Since the atomic velocities are nonrelativistic, the atomic momentum in Eq. (5) can be written classically p a = m 0 v a , where m 0 is the mass and v a is the velocity of the atom. The momentum density of the MDW, driven forward by light, is then given by 22 The energy density of the atomic MDW in the approximation of negligible kinetic energy of atoms is given by where n a0 is the atomic number density of the material before the optical force has started to displace the atoms.
In locally homogeneous regions of space, the momentum and energy densities of the atomic MDW, resulting from Eqs. (6) and (7) and the optical force density in Eq. (3), are given by 23 Here n p and n g are the local phase and group refractive indices of the material at the central angular frequency of the narrow frequency band. In Eq. (8), the angle bracket is given by This time-and position-dependent quantity is also equal to the time average of the EM momentum density over the harmonic cycle.
The total momentum density of the MP state of light is a sum of the EM momentum density in Eq. (1) and the momentum density of the MDW in Eq. (6) as Correspondingly, the total energy density of the coupled MP state of light is a sum of the EM energy density in Eq. (2) and the energy density of the MDW in Eq. (7) as See "Methods" for a brief comparison of the total momentum density of light in the MP theory in Eq. (10) with the total momentum density of light used by Veselago 10 .

Light in a hypothetical homogeneous NIM.
Here we apply the theory to homogeneous NIMs, even though such materials are not presently available. The total momentum and energy of the atomic MDW are volume integrals of the classical momentum and energy densities of atoms in Eqs. (6) and (7). For a homogeneous medium, one thus obtains 22 Here p EM = G EM d 3 r and E EM = W EM d 3 r are the total EM momentum and energy of the light pulse. The last forms of Eqs. (12) and (13) are straightforwardly obtained from Eqs. (8) and (9) or from the MP quasiparticle model for dispersive media 25 .
By summing the EM field and atomic MDW contributions, one obtains the total momentum and energy of the coupled MP state of light in homogeneous NIMs as www.nature.com/scientificreports/ The total MP momentum in Eq. (14) is directed opposite to the volume integral of the Poynting vector and the related total EM momentum for NIMs with n p < 0 , and it is discussed further in "Methods". One can conclude that the optical force density in Eq. (3) on average reduces the number density of atoms from n a0 to n a , and thus, makes the total energy of the coupled field-material state of light in Eq. (15) negative. Therefore, in the MP theory of light, the total energy and the component of the momentum in the direction of propagation of light are both negative in homogeneous NIMs.
We visualize the dynamics of the material and the emergence of the atomic MDW consisting of local rarefications of the material density by using a Gaussian light pulse and the simulation parameters of a hypothetical homogeneous NIM described in "Methods". The Gaussian light pulse has a central vacuum wavelength of 0 = 1746.7 nm, the peak electric field amplitude of E 0 = 2744.9 V/m, and the relative spectral width of � 0 / 0 = 0.015 . Figure 2a shows the optical force density of Eq. (3) for the simulated Gaussian light pulse at an instance of time when the center of the light pulse is at position x = 8 μm. Figure 2b presents the EM and atomic MDW momentum densities of Eqs. (1) and (6). Figure 2c depicts the EM and atomic MDW energy densities of Eqs. (2) and (7). The values of the EM quantities are positive at all positions and times while the MDW quantities are negative. The atoms are drawn backward and thus, the density of the material is locally rarefied. The atomic displacement r a = t −∞ v a dt ′ is depicted in Fig. 2d. It shows how the light pulse leaves atoms behind it displaced backward. The atomic displacement is very small, but it can be nanometers for pulses that are longer and have higher intensity as shown in Fig. 7 of Ref. 25 .
Consistency with the principles of relativistic mechanics. We next show by a quasiparticle analysis that the positive total energy density of light of earlier works 10,31,32 on homogeneous NIMs is very challenging since, in relativistic mechanics, the momentum per energy ratio must be given by p/E = v g /c 2 , where v g is the group velocity 11 . Using Eqs. (14) and (15) and the EM energy E EM = W EM d 3 r = ω 0 normalized for a singlephoton, the MP energy and momentum in a homogeneous NIM are given by 25 . Thus, the MP quasiparticle fulfills the standard momentum per energy ratio p/E = v g /c 2 of relativistic mechanics 11 . This ratio is violated in conventional theories of light in homogeneous NIMs 10,31,32 pinpointing the importance of the MDW energy included in the present theory.
Light in an inhomogeneous NIM. In inhomogeneous NIMs, i.e., in all NIM structures that have been experimentally realized so far, the total momentum and energy of light include momentum and energy terms related to the gradients of material parameters and internal material interfaces of the NIM structure. This means that the MDW in typical NIM geometries is split into the volume and internal interface contributions. Thus, for the MP state of light, we obtain These relations are more complex than the corresponding relations for a homogeneous NIM in Eqs. (14) and (15). Consequently, the values of p MP and E MP following from Eqs. (16) and (17) for general NIM structures can be either positive or negative. Whether it is possible to construct effective parameters in terms of which the total momentum and energy of light in an inhomogeneous NIM can be written using the EM quantities, as in Eqs. (14) and (15) for a homogeneous NIM, is an open question left as a topic of further work. In any case, it can be shown that the MP energy and momentum in Eqs. (16) and (17) fulfill the standard momentum per energy ratio p/E = v g /c 2 of relativistic mechanics 11 .
The internal interface MDW momentum and energy contributions p MDW in Eqs. (16) and (17) do not include contributions from the external interfaces of the NIM structure. The excluded external interface contributions p interface and E interface are left at the interface when the light pulse enters the NIM structure. They decay afterward by internal phonon relaxation of the structure. These terms also account for the partial reflection taking place at the interface. The terms p interface and E interface are essential in the conservation laws when light enters or exits the NIM structure. We also note that, in the present work, the conventional boundary conditions of the EM fields are used at all material interfaces. This assumes that the material boundaries are sharp. In general, this is not the case, and the surfaces can be rough, for example. If the structure size is made close to the atomic scale, the surfaces cannot be sharp, and it is necessary to account for the related surface effects.
An example of an inhomogeneous NIM structure is schematically illustrated in Fig. 3. The structure consists of dielectric cylinders embedded in a plasmonic host material. The propagation of light in this kind of a structure has been previously theoretically investigated, e.g., by Costa et al. 18 . In our setting, there are approximately seven unit cells in a single vacuum wavelength of the light pulse. More detailed material and structure parameters are shown in Fig. 3 and discussed in "Methods". The parameters of the incident Gaussian light pulse are the same as those in the case of a homogeneous NIM above. Figure 4 presents the optical volume and interface force densities resulting from the Gaussian light pulse. The snapshot is taken at the instance of time when the center of the Gaussian pulse is at the position of the entry interface of the NIM structure. Figures 5a-d show the position dependencies of the EM momentum density, (14) p MP =p EM + p MDW = n p n g p EM , www.nature.com/scientificreports/ www.nature.com/scientificreports/ MDW volume momentum density, interface momentum density, and the MDW mass density, respectively, at the instance of time corresponding to the force densities in Fig. 4. The EM momentum density in Fig. 5a consists of the incident and transmitted contributions. The reflection is close to zero for the NIM structure studied. The atomic MDW momentum density in Fig. 5b is zero on the left since there is no MDW in vacuum. It has the highest values in the dielectric cylinders, but it is also locally nonzero in the plasmonic material, where it is, however, so much smaller than in dielectric cylinders that this cannot be seen in the color scale of the figure. The interface momentum density in Fig. 5c is directed close to the normal of the interfaces. Part of this momentum density stays in the vicinity of the vacuum-NIM interface when the light pulse has passed while part of it propagates with light inside the NIM. Figure 5d presents the mass density of the atomic MDW, which is driven forward by the optical volume force density. Like the MDW momentum density in Fig. 5b, the MDW mass density is also locally small but nonzero in the plasmonic material. The MDW mass density follows the Gaussian envelope of the light pulse. Apart from different relative magnitudes in the dielectric cylinders and in the plasmonic host material, it reminds the EM energy density. The EM energy density is depicted in a larger scale in Fig. 5e at a later instance of time when the trailing edge of the pulse has also entered the NIM. At this instance of time, we have also plotted the spatially averaged EM, MDW volume, and internal interface momentum and energy densities, and the atomic displacement, in Fig. 5f-h, respectively. The spatial averaging is made over a distance of one wavelength. These figures allow for transparent comparison with the case of a homogeneous NIM in Fig. 2. It is seen that there are negative contributions in the momentum and energy densities, which originate from the internal interfaces in the MDW. However, in the present case, the volume contributions of the MDW are more positive and the net effect is a slight densification of the material, which contrasts to the rarefication seen in Fig. 2. From the negative internal interface contribution, it however, results that the net momentum propagating with the light pulse in the inhomogeneous NIM is smaller than the momentum of the pulse in vacuum. This is discussed in more detail for the integrals of the momentum density contributions below. Apart from the slight spreading of the pulse in the nonlinear dispersion of the inhomogeneous NIM, the spatially averaged EM momentum and energy densities in Fig. 5f,g are approximately equal to what would be obtained by correspondingly averaging the EM quantities of the homogeneous case in Fig. 2b,c. When the light pulse enters the NIM, the incident EM momentum component p 0 is split to partial momenta carried by the reflected and transmitted parts of the EM field, the atomic MDW, and the interface atoms. The comparison of these momenta between the homogeneous and inhomogeneous NIM cases studied above is presented in Table 1. The total MP momentum component of the coupled field-material state of light in both cases is smaller than p 0 . In the homogeneous case, it is negative while, in the inhomogeneous case, it is positive. Since the total MP momentum component is smaller than p 0 , the interface momentum is directed toward the NIM due to the conservation law of momentum. The same applies at the interface where the light pulse exits the NIM. Experimental detection of reduced or negative radiation pressure. A schematic experimental setup for probing radiation pressure of light in NIMs is illustrated in Fig. 6a. The setup enables the measurement of the optical force produced when light exits the NIM structure. The atomic MDW densification or rarefication waves are also present in this setup which, however, focuses on the force measurement [33][34][35][36][37][38][39][40] . In contrast to the NIM in Fig. 3, here we use a NIM where the host material is vacuum and there are cylinders and c-shaped rods made of a plasmonic material. See Refs. 4,9 for examples of artificial structures that can be fabricated with presently available technologies. The NIM is split into bulk-NIM and layer-NIM parts. The layer-NIM is attached to a cantilever, which is able to bend. Therefore, the force experienced by the layer-NIM can be measured by detecting the bending of the cantilever under the continuous-wave illumination of the NIM. This bending can be detected with a probe laser beam and a photodetector utilizing the optics of an atomic force microscope. Cantilever displacements of fractions of a nanometer are measurable. The setup of Fig. 6a is not optimized for www.nature.com/scientificreports/ the actual measurement but it is meant to provide an idea of the measurement that can be realized with various kinds of NIM structures. Figure 6b presents the time-averaged EM momentum density. It shows that, on average, the energy and momentum fluxes of the EM field are directed forward in the cylinders and, counterintuitively, backward in the c-shaped rods. Figure 6c depicts the force density averaged both over the harmonic cycle and over vertical layers of the NIM structure. The layer-NIM attached to the cantilever experiences the strongest force and it is directed backward, i.e., the radiation pressure is reduced or negative. Qualitatively, the same result applies to various NIM Table 1. Momentum components after the light pulse has entered the NIM. When the light pulse enters the NIM, the incident EM momentum component p 0 in the direction of propagation is split to partial momentum components carried by the reflected EM field ( p  www.nature.com/scientificreports/ structures in which the total momentum of light is smaller than that of vacuum, i.e., the momentum of light does not need to be negative. The force directed inward the NIM contrasts the force in normal materials where it is directed outward 41 . Since the host material of the NIM used in the setup is vacuum, the conventional Maxwell stress tensor of vacuum could also be applied giving results equal to those of us.

Conclusions
In conclusion, in contrast to normal materials, the total momentum of light in NIMs can be smaller than the momentum of light in vacuum. The radiation pressure in NIMs can even become negative depending on the subwavelength structure of the material. In NIMs exhibiting negative radiation pressure, the negative total momentum of light results from the atomic mass density rarefication waves caused by the optical force density. www.nature.com/scientificreports/ In the corresponding single light quantum picture, the MP quasiparticles have a negative total energy, which allows the total momentum of light to be directed opposite to the group velocity. The backward atomic velocities and the resulting mass density rarefication are real, experimentally measurable, and necessary for the fulfillment of the constant center-of-energy velocity of an isolated system. The theory applies to both hypothetical homogeneous and conventional inhomogeneous NIMs structures. We have presented a schematic experimental setup, which allows the measurement of the reduced or negative radiation pressure effect of light in NIMs. The present results demonstrate the versatile use of the MP theory of light for modeling optomechanics of light in general photonic crystal structures.

Methods
Energy densities. In a lossless material, the energy density of the EM field can be solved from the local conservation law of energy, 1 where the right-hand side can be approximated by zero due to its smallness, giving W EM = −c 2 t −∞ ∇ · G EM dt ′22,42 . This integral expression then leads to the result, given in Eq. (2) 23 . Solving the energy density of the MDW from the conservation law of energy correspondingly gives Eq. (7), which for homogeneous regions of space is equal to Eq. (9) 23 .
Total momentum density of light. The instantaneous momentum density of the MP state of light in Eq. (10) is not equal to the total momentum density of light used by Veselago 10 , i.e., n 2 p E × H/c 2 + 1 2 k( ∂ε ∂ω 0 E 2 + ∂µ ∂ω 0 H 2 ) , where k is the wave vector of light in the material. For homogeneous NIMs, one can, however, show that the harmonic cycle time averages of these momentum densities are equal. Veselago did not consider the possibility that part of the total momentum of light could be carried by the material. Therefore, his approach deviates from the MP theory. For homogeneous NIMs with n p < 0 , the MP momentum in Eq. (14) is directed opposite The NIM crystal is fabricated on a substrate so that the cylinders and c-shaped rods are fixed from their ends (not shown). The NIM crystal is then split into bulk-NIM and layer-NIM parts. The bulk-NIM is attached to a rigid supporting structure while the layer-NIM is attached to a cantilever. The gap between the bulk-and layer-NIMs is exaggerated to make it visible. In the proposed experiment, this gap is equal to the gaps inside layers of the bulk-NIM. The cantilever is able to bend as a result of the optical force. The nanoscale bending ( ≪ a ) of the cantilever is recorded by monitoring the reflection of a probe laser beam using a photodetector. The right inset shows the NIM structure used in our example. (b) Time-averaged EM momentum density of a continuous wave light beam propagating through the NIM. The linearly scaled color indicates the magnitude, and the logarithmically scaled arrows indicate the direction. (c) The x component of the total optical force density averaged over the harmonic cycle and over the vertical layers of the NIM. The bulk-NIM is at x < 0 and the layer-NIM is at x > 0. www.nature.com/scientificreports/ to the EM energy and momentum flux. This fundamentally originates from the backward-directed optical force in Eq. (3) and the resulting backward-directed atomic velocities and the MDW momentum density. For normal dispersive materials with n p > 0 , the total momentum of the MP state of light in Eq. (14) explains the results of the high-precision measurements of radiation pressure by Jones et al. 43 .

MP quasiparticles.
In the single light quantum picture, the energy and momentum of the MP quasiparticle can be shown to be equal to the volume integrals of the classical energy and momentum densities normalized for the known single-photon EM energy E EM = W EM d 3 r = ω 0 . For homogeneous NIMs, this leads to E MP = W MP d 3 r = n p n g ω 0 and p MP = G MP d 3 r = (n p ω 0 /c)v g /|v g | 22,25 . Using the relativistic energymomentum relation, given by (m MP c 2 ) 2 = E 2 MP − (p MP c) 2 , one then obtains the MP quasiparticle rest mass as m MP = n p n 2 g − 1 ω 0 /c 2 . Therefore, for homogeneous NIMs, m MP is negative. This result can also be obtained by using the relativistic MP quasiparticle model presented in Ref. 25 . The energy and momentum of the MP quasiparticle in all inertial frames are given by the same relativistic relations as for all other particles and quasiparticles with a nonzero rest mass. Thus, for the MP state of light, we have E MP = γ v g m MP c 2 and p MP = γ v g m MP v g , where v g is the group velocity and γ v g = 1/ 1 − |v g | 2 /c 2 is the corresponding Lorentz factor.
Gaussian light pulse. The electric field of a one-dimensional Gaussian light pulse 44 , linearly polarized in the y direction and propagating in the x direction in vacuum, is given by E(r, t) = E 0 cos[k 0 (x − ct)]e −(�k 0 ) 2 (x−ct) 2 /2ŷ . Here ŷ is the unit vector in the direction of the y axis, E 0 is the electric field amplitude, k 0 = ω 0 /c = 2π/ 0 is the wave number, and k 0 is the standard deviation of the wave number. In all simulations of the present work, we use E 0 = 2744.9 V/m, k 0 = 0.015k 0 , and 0 = 1746.7 nm. The wavelength is chosen so that it corresponds to the case when the effective phase refractive index of a realistic NIM studied in Ref. 6 becomes n p,eff = −1 . We have selected the NIM structures of the present work so that this condition is also satisfied for them. The magnetic field corresponding to the electric field above is determined by Maxwell's equations. The fields inside materials also follow from Maxwell's equations when the pulse penetrates through the vacuum-material interface. The simulations are carried out using Comsol Multiphysics simulation tool 45 . Periodic boundary conditions are applied in the direction of the y axis in Fig. 3.

Dispersion relation of the EM field.
In the illustration of light in a homogeneous NIM, we use a dispersion relation linearized around the central angular frequency as ω(k) = ω 0 + (c/n g )(k − n p k 0 ) . This dispersion relation corresponds to the frequency-dependent phase refractive index n p (ω) = n g + (n p − n g )ω 0 /ω . The group refractive index is constant as n g (ω) = ∂ ∂ω [ωn p (ω)] = n g . This means that the light pulse envelope preserves its shape while the pulse is propagating. To make reflection of light from the vacuum-NIM interface negligible, we match the wave impedance of NIM to that of free space byassuming equal frequency dependencies for the permittivity and permeability by writing ε(ω) = ε 0 n p (ω) and µ(ω) = µ 0 n p (ω) . These are both negative at ω 0 if n p is negative. Thus, the condition for negative refraction, i.e., ε|µ| + µ|ε| < 0 46 , is satisfied in our example. For the phase and group refractive indices at 0 = 1746.7 nm, we use the values of n p = −1 and n g = 10.52 corresponding to the central frequency values of the effective refractive indices calculated for the inhomogeneous NIM. The plasmonic host material of the inhomogeneous NIM in Fig. 3 is assumed to obey the Drude model dispersion relation with a relative permittivity equal to ε h /ε 0 = 1 − ω 2 p ω(ω+iŴ) . Here the plasma frequency is ω p = 1.22573 × 10 15 rad/s. The collision frequency Ŵ is set to zero so that the EM energy is conserved and we can study the field-driven MDW, which is separate from the impulse taken by atoms in absorption. The dielectric cylinders in Fig. 3 have a constant relative permittivity of ε c /ε 0 = 50.47 , which is the value used for a similar structure in Ref. 18 . The plasmonic material used in the schematic experimental setup of Fig. 6 has a plasma frequency of ω p = 1.2459 × 10 15 rad/s.

Mechanical constants of the materials.
In the simulations, we have used the equilibrium mass density of ρ a0 = 2200 kg/m 3 for the solid materials. The elastic force density was neglected in the present time-dependent simulations since its effect is very small in the short time scale of the propagation of the light pulse 24,25,47 . The elastic force density can be added in Newton's equation of motion without any change in the theory of the coupled field-material dynamics. It is important if one wants to analyze the relaxation of the nonequilibrium atomic displacements produced by the optical force density of the light pulse. This relaxation takes place in the longer time scale of sound waves. In the case of continuous-wave illumination in the experimental setup in Fig. 6, the role of elasticity is to balance the optical force on the layer-NIM by the opposite reaction force resulting from the bending of the cantilever. Thus, in this example, the elastic force density should be added in the simulations if one wants to simulate the bending of the cantilever. In this work, we restrict our simulations of this setup to the calculation of the optical force.

Data availability
The simulation results presented in this paper are available from the corresponding author on reasonable request.